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Abstract. In the non viscous fluid dynamics, Smooth Particle Hydrody- 
namics (SPH), as a free Lagrangian "shock capturing" method adopts either an 
artificial viscosity contribution or an appropriate Riemann solver technique. An 
explicit or an implicit dissipation, introduced in such techniques, is necessary to 
solve the Euler equations to solve flow discontinuities (the Riemann problem). 
Dissipation is useful to smooth out spurious heating and to treat transport phe- 
nomena. A simple, effective solution of the Riemann problem is here proposed, 
based on an empirical reformulation of the equation of state (EoS) in the Euler 
equations in fluid dynamics, whose limit for a motionless gas coincides with the 
classic EoS of ideal gases. Results on ID shock tube tests are here shown, as 
well as a 3D transport application on accretion discs in close binaries (CBs). 



1. Introduction 



A dissipation is normally adopted to handle discontinuities in non viscous flows 
(the Riemann problem). Such a dissipation can be either explicit (artificial vis- 
cosity) or implicit, being intrinsic within the adopted algorithm (e.g. Godunov- 
like methods). An artificial viscosity is introduced in SPH, as a shock capturing 
method, to prevent particle interpenetration and to smooth out spurious heat- 
ing in the solution of strictly hyperbolic system of Euler equations. Such a 
small dissipation also produces both mass and angular momentum transport 



in SPH physically inviscid modelling of accretion disc s (Moltcni et al. 



1994; Lanzafame et al. 



Chak rabarti fc Molteni 
19971 



1991 119931. 11991 l200d. 12001]: iBelvedere et al. 



1991 



1993; 



1993; Meglic ki et al. Hl993l : lMurrav lll996l : lLanzafame fe Belvedere 



Morris fc Monaghan 



j, .19981 : lOkazaki et al. 1 12002 ) . Efforts were accomplished in SPH to solve 
both the " approximate" and the " exact" Riem ann problem, ei t her via an explicit 
reformulation of the artific ial viscosity term (Bal saral Il995l ; iMonaghan" 1997; 



199 7]) or via sophisticated G o dunov algorithms ( Yukawa et al. 



19971 : llnutsuka feooaTPajshikov fe Medin I [2001 : iMolteni k Bilello 1 120031 ) . A 



physical way to solve the Riemann problem, based on a reformulation of the EoS 
in the Euler equations, is here presented, where particle SPH pressure terms are 
recalculated without any artificial viscosity contribution. Since shock flows are 
non equilibrium events, we pay attention to the fact that the EoS: p = (7 — l)pe 
for ideal flows cannot exactly be applied to solve fl ow discont inuities. Applica- 
tions of the reformulated EoS to ID shock tubes dSod | [T978). to solve shocks 



are shown, as well as a comparison to models adopting a different dissipation is 
given for the coming out of spiral patterns in accretion discs in CBs. 
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2. The Euler equations: SPH formulation 

In the Lagrangian ideal non-viscous gas hydrodynamics, the relevant equations 
(Euler equations) are: 

+ pV ■ v = continuity equation (1) 

dt 

dv Vp 

—— = momentum equation (2) 

dt p H w 

de p 

— = V • v energy equation (3) 

dt p 

p = (7 — l)pe perfect gas equation (4) 
dv 

-= = v kinematic equation (5) 

The most of the adopted symbols have the usual meaning: d/dt stands for 
the Lagrangian derivative, p is the gas density, e is the thermal energy per unit 
mass. The adiabatic index 7 has the meaning of a numerical parameter whose 
value lies in the range between 1 and 5/3, in principle. 

The SPH method is a free Lagrangian scheme that discretizes the fluid into 
moving interacting and interpolating domains called "particles". All particles 
move according to pressure and body forces. The method makes use of a Kernel 
W useful to interpolate a physical quan t ity A (r) related to a gas particle at 
position r according to (IMonaghan I ll98.1ll992h : 



r N A ■ N A- 

A{r) = / A(r')W(r,r', h)dr' ^A i = J2 —W^rj, h) = J2 ~ w ij ( 6 ) 
Jd - =l rij . =1 rij 

W(r,r', h), the interpolation Kernel, is a continuous function - or two con- 
necting continuous functions whose derivatives are continuous even at the con- 
necting point - defined in the spatial range 2h, whose limit for h — > is the 
Dirac delta distribution function. All physical quantities are described as ex- 
tensive properties smoothly distributed in space and computed by interpolation 
at r. In SPH terms the sum is extended to all particles included within the 
interpolation domain D, rij = Pj/rrij is the number density relative to the jth 
particle. Wir^^ph) < 1 is the adopted interpolation Kernel whose value is 
determined by the relative distance between particles i and j. 

In SPH formalism, equations (2) and (3) take the form, respectively: 



dUi _ \- p\ , Pj 



dt 

N 



2 




ViWij (7) 



where = — Vj and rrij is the mass of jth particle. 
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For a better energy conservation, the total energy E = (e + \v 2 ) can also 
be introduced in the SPH formulation: 



^ = -E^ + ^ 'V# y (9) 



= - E m i + —> 
dt V pi Pj j 

of the energy equation: d{e + v 2 /2)/dt = — V • (pv)/p. 
In this scheme the continuity equation takes the form: 

, N 

^ = E™^- V ^ (io) 

or, as we adopt, it can be written as: pi = J2jLi m jWij, which identifies 
the natural space interpolation of particle densities according to equation (6). 

The pr e ssure terms p* i nclude the artificial viscosity co ntribution given by 
Monaehan I (Il985l . I1992T ) and M onaghan fe Lattanzio I (j 19851 ) , with an appropri 



ate thermal diffusion term which reduces shock fluctuations: p* = pi(l + rjij). 
r]ij is given by: 



+ 0p 2 j , where 



2hv . . -r ■ ■ 
— »j — »j 



otherwise 



(11) 



with c s i being the sound speed of the ith particle, r 



■'j 



T | 



■r 



e « h 2 , 

a ~ 1 an d 2. These a and parameters of the order of the unity are usually 
adopted ( Lattanzio et al. Ill985l ) to damp oscillatio ns past high Mach n umber 
shock fronts developed by non- linear instabilities ( Boris Sz Book 19731 ). The 
linear a term is based on the viscosity of a gas. The quadratic (0, Von Neumann- 
Richtmyer-like) artificial viscosity term is necessary to handle strong shocks. 
Linear a and quadratic artificial viscosity terms are ~ 1. A reformulation 
of the artificial viscosity could be even necessary because, for "weak shocks" or 
low Mach numbers, the fluid becomes "too vi scous" and ang u lar momentum and 
vorticity could be non-physically transferred. Bal sara ( 19951 ) : iMonaghan (1997) 
developed some effective techniques to get a limitatio n of the effectiveness of the 
artificial viscosity, while iMorris Sz Monaghan" ( 19971 ) derived a new formulation 
for artificial viscosity conceptually based on the Riemann problem. 

The solution of the Riemann problem in SPH is also obtained at the inter- 
particle contact points among particles, where a pressure and a velocit y, relative 



to the f low discontinuity, are comp ut ed. This is also cle a rly sh own in llnutsuka 



(2002); iParshikov k Medin I (120021 '): iMolteni fc Bilello I (jl003), where the new 
pressure p* and velocity v* are reintroduced in the Euler equations to obtain 
the new solutions compatib le with inyiscid flow discontinuities. In S PH, we 
pay attention in particular ( Inutsuka 2002 : Parshik ov Sz Medin 20021 ) to the 
particle pressures p* and p*j, in the SPH formulation of the momentum and en- 
ergy equations (7) and (8 or 9), whose substitution with pressures, solutions 
of the Riemann problem, excludes any artificia l viscosity adoption, although a 
dissipation is necessarily implicitly introduced ( Park Sz Kwon [20031 ). 
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3. How EoS matches the Riemann problem 

The classical EoS for a perfect gas: p = (7 — l)pe is strictly applied in fluid 
dynamics when the gas components do not collide with each other whenever a 
thermal equilibrium exists. It modifies in: p* = (7 — l)pe+ other in the case of 
gas collisions . The further ter m takes into account the velocity of perturbation 
propagation ( Monaghan II 19971 ). This velocity equals the ideal gas sound velocity 



c s whenever we treat static gases or in the case of rarefaction waves. Instead, 
it includes the "compressio n velocity": •£ ,;,-/ | r_ij I in the case of shocks 
(v s ig^j = c s i — v_i ■ tjj 1 1 Ljj I ; I Whitehurst I (Il995l )). In the first case, we write the 



EoS for inviscid ideal gases as: p = pc 2 /^, where c s = ( r yp/p) 1 ^ 2 = [7(7 — l)^] 1 ^ 2 - 
In the second case, the new formulation for the EoS is obtained squaring v s i g ^j, 
so that c 2 s {l — v s hock/cs) 2 is an energy per unit mass in the case of shock. Hence: 

£ c 2(l_l^) 2 if Vshock <Q 
,>-[ V " (12) 

^c 2 if v shock > 



In the SPH scheme, being: 



7 " 



if V;a ■ Taa < 



* Pi 2 (a v shock,i\ I fcjl V V 

Pi = ~ c si 1 > V s hock,i = \ (13) 



otherwise. 



This formulation introduces the "shock pressure term" p{v 2 hock — 2w s fc oc fcC s )/7, 
whose dependence on linear and quadratic power on ■ r^/ \ r_ij I is analogue 
to both the linear and the quadratic components of the artificial viscosity term 
(eq. 11). The linear term in -r^ / j \ is based on the viscosity of a gas. The 
quadratic term in v_ij • Lij / \ Lij \ (Von Neumann - Richtmyer - like) is necessary 
to handle strong shocks. These contributions involve a dissipative power, whose 
effect correspond to an increase of the gas pressure. Therefore, we adopt the 
formulation (eq. 13) as p* in the SPH formulation of the momentum (eq. 8) 
and energy equations (eqs. 9 or 10): 

ViWy (14) 
•ViWy, (15) 




ViWij. (16) 

However, to give a generalization, we need only one general EoS and not a 
separation of the EoS according to the kinematic of the flow. To this purpose, 
we can generalize the EoS: p* = pc 2 {l — v s h oc k/c s ) 2 h as: 
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Figure 1. ID shock tube tests as far as both analytical (solid line) and 
our SPH-Ricmann (dots) results are concerned. Density p, thermal energy e 
and velocity v are plotted at time T — 100. Density and thermal energy of 
particles initially at rest at time T — are also reported. The initial velocity 
is zero throughout. 



where C — > 1 for vr = ■ r^/ \ r_%j |< 0, whilst C — > otherwise. A 
simple empirical formulation can be: C = &rccot(RvR/ c s ) /ir , where R S> 1. R 
is a large number describing how much the flow description corresponds to that 
of an ideal gas. To this purpose, R w X/d, being A oc /? -1 / 3 the molecular mean 
free path, and d the mean linear dimension of gas molecules. 



4. Tests: ID Sod shock tube and 3D accretion discs in CBs 

A comparison of analytical and our SPH - Riemann ID shock tube test results 
(Sod 1978), also with the initial particle configuration (time T = 0), is made. 
Figg. 1 and 2 shows results concerning the particle density, thermal energy 
per unit mass and velocity at the same final computational time (T = 100). 
Throughout our SPH- Riemann simulations, the initial particle resolution length 
is h = 5 ■ 10~ 2 . The whole computational domain is built up with 2001 particles 
from X = to X = 100, whose mass is different, according to the shock initial 
position. At time T = all particles are motionless. 7 = 5/3, while the ratios 
P1/P2 = 3 and e\/e2 = 2 (Fig. 1 - left side), and P1/P2 = 3 and €\/e2 = 1 
(Fig. 1 - right side), between the two sides left-right. The first and the last 
5 particles of the ID computational domain, keep fixed positions and do not 
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Figure 2. XY plots of 64 greytones density p isocontours of 3D disc mod- 
elling in CBs. SP H refers to norm al SPH results. Disc model "a" refers to SPH 
modelling where LBalsara J (|1995t) attenuatio n is introduced. Model "b" refers 
to SPH modelling where iMonaghan I (|1997l ) treatment of the Riemann solver 
i s intr oduced. Model "c" refers to SPH modelling where [Morris fc Monaghanl 
(|1997f ) a decay is introduced. Model "d" refers to SPH modelling where the 
EoS reformulation is introduced. 



move. Discrepancies involve only 4 particles at most. This means that the 
physical dissipation introduced in the EoS (eq. 13) is effective. 

Results on 3D SPH simulations of accretion disc structure and dynamics in a 
CB are shown as far as the coming out of spiral patterns and shocks in the radial 
flow are concerned. The resolution length of SPH particles is h = 0.005, being 1 
the non-dimensional separation d±2 = 10 6 Km of the two stars. The injection of 
particles from the inner Lagrangian point LI is supersonic: Vi n j = 130 Km s , 
whilst the temperature of gas coming from the secondary star is T = 10 4 K and 
7 = 1.01. The compact primary is a 1M Q star, while the donor companion is 
a 0.5Mq star. The comparison with other different SPH techniques (see Fig. 
2 caption), where artificial dissipation is explicitly introduced, shows that the 
coming out of spirals, as well as other disc geometric details, are much better 
evidenced in those simulations where the EoS and related dissipation are treated 
in their full physical sense. 

3D accretion disc simulations clearly show different results. The adoption 
of a supersonic injection kinematics from LI is justified in terms both of con- 
servation of the flux momentum and in terms of t he Jacobi constant and in 
terms of the Bernoulli's theorem ( Lanzafame 1 1200*91 ) . This kinematics concern 
active phases of accretion discs where both radial disc increases and spiral struc- 
tures develop. Causes generating spiral patterns and disc geometric details are 
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deeply discussed in Sawada et al. (19871"); Smv 


it et al. (1987): Kaisie ( 


1989): 


Sawada Sz Matsuda ( 19921); ISavonniie et al. ( 


1994): Lanzafame et al. (2000. 


200lh as a consequence of external perturbations as tidal torques as well as the 



collisional impact of the injected flow onto the disc's outer edge. Hence, the 3D 
simulated structure better corresponding to such theoretical and observational 
constraints is that produced according to the formalism here proposed (model 
d). The typical 3D SPH disc structure (model a) does not fully correspond to 
these characteristics, while other 3D disc models lack of several features. 
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